/*
 * System.cpp
 *
 * Created on: 23 May 2011
 *   Author: allan
 */

#include "System.h"

// C++ includes
#include <algorithm>
#include <cmath>
using namespace std;

// Boost includes
#include <boost/foreach.hpp>

System::System()
{}

System::System(const vector<PhasePtr>& phases, const vector<ReactionPtr>& reactions) :
MultiPhase(phases), reactions(reactions)
{
	set<String> primarySpeciesSet, secondarySpeciesSet;
	
	// Loop over all the reactions
	BOOST_FOREACH(const ReactionPtr& reaction, reactions)
	{
		secondarySpeciesSet.insert(reaction->GetReactant().name);
		
		BOOST_FOREACH(const Reaction::Species& product, reaction->GetProducts())
			primarySpeciesSet.insert(product.name);
	}
		
	primarySpecies.assign(primarySpeciesSet.begin(), primarySpeciesSet.end());
	secondarySpecies.assign(secondarySpeciesSet.begin(), secondarySpeciesSet.end());
}

/// --- Getter methods --- ///
const StringVector& System::GetPrimarySpecies() const
{
	return primarySpecies;
}

const StringVector& System::GetSecondarySpecies() const
{
	return secondarySpecies;
}

const vector<ReactionPtr>& System::GetReactions() const
{
	return reactions;
}

//	/// --- Add methods --- ///
//	void
//	Open();
//	
//	void
//	AddPhase(const GenericPhase& phase);
//	
//	void
//	AddReaction(const GenericReaction& reaction);
//	
//	void
//	Close();

//void System::Open()
//{
//	// Open the phase assemblage object to allow the insertion of phases
//	multiphase.Open();
//}
//	
//void System::AddPhase(const GenericPhase& phase)
//{
//	phases.push_back( PhasePtr(phase.Clone()) );
//}
//
//void System::AddReaction(const GenericReaction& reaction)
//{
//	reactions.push_back( PhasePtr(reaction.Clone()) );
//}
//
//void System::Close()
//{
//	// Close the phase assemblage object
//	multiphase.Close();
//	
//	// Process all the inserted phases and reactions
//	
//	// The names of all species, which has been assembled from all the phases by the phase assemblage object 
//	const StringVector speciesNames = multiphase.GetSpecies();
//	
//	// A SpecieMap object to retrieve species index; used to set up the indexes data of the reactions
//	const SpeciesMap speciesMap(speciesNames);
//	
//	// Loop over all the reactions and set the indexes of their species
//	BOOST_FOREACH(const ReactionPtr& reaction, reactions)
//		reaction->SetSpeciesIndexes(speciesMap);
//	
//	// Loop over all the species (names of them) to create the Species objects
//	uint index = 0; BOOST_FOREACH(const String speciesName, speciesNames)
//	{
//		// Find, among all provided (reactions), the one that has the reactant named (speciesName)
//		const ReactionPtr reaction =	
//				*std::find_if(reactions.begin(), reactions.end(), [&](const ReactionPtr& r) { return r->GetReactant().name == speciesName; } );
//		
//		// Create a System::Species object
//		const Species newSpecies = { speciesName, index, reaction };
//		
//		// Store the new species
//		species.push_back(newSpecies);
//		
//		// Update the index variable
//		++index;
//	}
//}
//
//System::System()
//{}
//
//System::System(const vector<PhasePtr>& phases, const vector<ReactionPtr>& reactions) : 
//MultiPhase(phases), reactions(reactions)
//{
//	// Set the primary species of the system
//	set<Index> primarySpeciesSet;
//	
//	// Loop over all the reactions in the system to set up the indexes of the primary species, which are the products of the reactions
//	BOOST_FOREACH(const ReactionPtr reaction, reactions)
//		BOOST_FOREACH(const Reaction::Species& product, reaction->GetProducts())
//			primarySpeciesSet.insert(product.index);
//	
//	primarySpecies.assign(primarySpeciesSet.begin(), primarySpeciesSet.end());
//	
//	// Loop over all the species in the system to set up the indexes of the secondary species
//	for(uint ispecies = 0; ispecies < species.size(); ++ispecies)
//		if(primarySpeciesSet.find(ispecies) == primarySpeciesSet.end()) // Check if the current species index is NOT in the set of primary species
//			secondarySpecies.push_back(ispecies);
//	
//	// Setup the reaction matrix between secondary and primary species
//	vij = StoichiometricMatrix(secondarySpecies);
//}
//
//const System::ReactionRange System::GetReactions() const
//{
//	auto begin = boost::make_indirect_iterator(reactions.begin());
//	auto end   = boost::make_indirect_iterator(reactions.end());
//	
//	return boost::make_iterator_range(begin, end);
//}
//
//const StringVector& System::GetPrimarySpecies() const
//{
//	return primarySpecies;
//}
//
//const StringVector& System::GetSecondarySpecies() const
//{
//	return secondarySpecies;
//}
//
//const MatrixXd System::StoichiometricMatrix(const IndexVector& reactionIndexes) const
//{
//	// The number of reactions (Nr) and primary species (Nj)
//	const uint Nr = reactionIndexes.size();
//	const uint Nj = primarySpecies.size();
//	
//	// The stoichiometric matrix for the set of reactions
//	MatrixXd vrj(Nr, Nj);
//	
//	// Loop over all the reactions
//	for(uint r = 0; r < Nr; ++r)
//	{
//		// Loop over all the primary species
//		for(uint j = 0; j < Nj; ++j)
//		{
//			const vector<GenericReaction::Species>& products = reactions[reactionIndexes[r]].GetProducts();
//			
//			// Find the (local) index of the product that corresponds to the current primary species j (note that the current reaction may not have this primary species as a product)
//			const uint iproduct = std::find_if(products.begin(), products.end(),
//					                  [&](const GenericReaction::Species& product) { return product.index == primarySpecies[j]; }) - products.begin();
//			
//			// Check if iproduct is within the limits [0, numProducts)
//			vrj(r, j) = (iproduct < products.size()) ? products[iproduct].stoichiometry : 0.0;
//		}
//	}
//	
//	return vrj;
//}
//
//void PrintWithIndexes(ostream& out, const StringVector& vec)
//{
//	uint index = 0; 
//	BOOST_FOREACH(const String s, vec)
//	{
//		out << "(" << index << ") " << s << ", ";
//		++index;
//	}
//}
//
//void PrintWithIndexes(ostream& out, const IndexVector& speciesIndexes, const System& system)
//{
//	BOOST_FOREACH(Index index, speciesIndexes)
//		out << "(" << index << ") " << system.GetSpecies()[index] << ", ";
//}
//
//ostream& operator<<(ostream& out, const System& system)
//{
//	uint iPhase = 0; 
//	BOOST_FOREACH(const GenericPhase& phase, system.GetPhases())
//	{
//		out << "Phase " << iPhase << ": " << phase.GetPhaseName() << endl;
//		out << "  -PhaseSpecies: " << phase.GetSpecies() << endl << endl;
//		
//		++iPhase;
//	}
//	
//	out << "Species:" << endl; 
//	out << "  -GlobalSpecies: " << endl; 
//	out << "      "; PrintWithIndexes(out, system.GetSpecies()); out << endl << endl;
//	
//	out << "  -PrimarySpecies: " << endl; 
//	out << "      "; PrintWithIndexes(out, system[system.GetPrimarySpecies()], system); out << endl << endl;
//	
//	out << "  -SecondarySpecies: " << endl; 
//	out << "      "; PrintWithIndexes(out, system[system.GetSecondarySpecies()], system); out << endl << endl;
//	
//	uint iReaction = 0;
//	BOOST_FOREACH(const GenericReaction& reaction, system.GetReactions())
//	{
//		out << "Reaction " << iReaction << ": " << endl; 
//		out << "   " << reaction << endl << endl; 
//		
//		++iReaction;
//	}
//	
//	out << "Stoichiometric Matrix vij: " << endl;
//	out << system.vij << endl << endl;
//	
//	return out;
//}
